/*
http://www.gnu.org/software/gsl/manual/html_node/GSL-BLAS-Interface.html -->Epe3igisi
http://www.gnu.org/software/gsl/manual/html_node/GSL-CBLAS-Library.html	--> antistixisi se CBLAS
*/

while( cblas_dnrm2 (n, r, 1)/cblas_dnrm2 (n, b, 1) > ITOL && iter < h)
	++iter;

	//SOLVE Mr = z;	-> r = inv(M)*z
	
	//theloume na diatirisoume ta dedomena tou r
	//den mas endiaferoun ta dedomena tou z prin tin pra3i
	//default pra3i me 1 stin diagonio
	cblas_dcopy (n, r, 1, z, 1);
	cblas_dtrsv (CblasRowMajor, CblasUpper, CblasNoTrans,CblasNonUnit, n, MNA[0], n, z, 1);

	//rho = dotproduct(r,z)
	double rho = cblas_ddot (n, r, 1, z, 1);
	
	if (iter == 1){
		//p = z
		cblas_dcopy (n, z, 1, p, 1);
	}
	else{		
		//p = z + beta*p;
		double beta = rho/rho1;
		
		cblas_dscal (n, beta, p, 1);//rescale the vector p by the multiplicative factor beta	
		
		cblas_daxpy (n, 1, z, 1, p, 1);//p = 1*x + p	
	}
	double rho1 = rho;
	
	//q = Aq
	cblas_dcopy (n, p, 1, q, 1);
	cblas_dgemv (CblasRowMajor,CblasNoTrans,n,n,1,MNA[0],n,q,1,0,q,1);	//q = 1*MNAxq + 0*q

	//alpha = rho/dotproduct(p,q)
	double aplha = rho/cblas_ddot (n, p, 1, q, 1);
	
	//x = x + aplha*p;
	cblas_dcopy (n, p, 1, temp_p, 1);
	cblas_dscal (n, alpha, temp_p, 1);//rescale the vector p by the multiplicative factor aplha
	cblas_daxpy (n, 1, temp_p, 1, x, 1);// sum x = 1*x + temp_p

	//R = R - aplha*Q;
	cblas_dcopy (n, q, 1, temp_q, 1);
	cblas_dscal (n, -alpha, temp_q, 1);//rescale the vector p by the multiplicative factor -aplha
	cblas_daxpy (n, 1, temp_q, 1, r, 1);// sum r = 1*r - temp_p











}
